{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multivariable Stability analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we used proportional control on SISO systems we observed that there is usually an upper bound on the controller gain $K_c$ above which the controlled system becomes unstable. Let's investigate the equivalent calculation for MIMO systems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy\n",
    "sympy.init_printing()\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = sympy.Symbol('s')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This matrix is from example 16.2 in Seborg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}\\frac{2}{10 s + 1} & \\frac{3}{2 \\left(s + 1\\right)}\\\\\\frac{3}{2 \\left(s + 1\\right)} & \\frac{2}{10 s + 1}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡   2           3    ⎤\n",
       "⎢────────   ─────────⎥\n",
       "⎢10⋅s + 1   2⋅(s + 1)⎥\n",
       "⎢                    ⎥\n",
       "⎢    3         2     ⎥\n",
       "⎢─────────  ──────── ⎥\n",
       "⎣2⋅(s + 1)  10⋅s + 1 ⎦"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Gp = sympy.Matrix([[2/(10*s + 1), sympy.Rational('1.5')/(s + 1)],\n",
    "                   [sympy.Rational('1.5')/(s + 1), 2/(10*s + 1)]])\n",
    "Gp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "K_c1, K_c2 = sympy.symbols('K_c1, K_c2', real=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unlike in SISO systems, we now have a choice of pairing. We will see that there are differences in the stability behaviour for the different pairings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "diagonal = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "if diagonal:\n",
    "    Gc = sympy.Matrix([[K_c1, 0],\n",
    "                       [0, K_c2]])\n",
    "else:\n",
    "    Gc = sympy.Matrix([[0, K_c2],\n",
    "                       [K_c1, 0]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "I = sympy.Matrix([[1, 0],\n",
    "                  [0, 1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The characteristic equation can be obtained from the $|I + GpGc|$. I divide by 4 here to obtain a final constant of 1 like in the example to make comparison easier. Make sure you understand that any constant multiple of the characteristic equation will have the same poles and zeros."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "charpoly = sympy.poly(sympy.numer((I + Gp*Gc).det().cancel())/4, s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare with Equation 16-20:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "charpoly2 = sympy.poly(\n",
    "    sympy.numer(\n",
    "        ((1 + Gc[0,0]*Gp[0,0])*(1 + Gc[1,1]*Gp[1,1]) - Gc[0,0]*Gc[1,1]*Gp[0,1]*Gp[1,0]).cancel()\n",
    "    )/4, s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "charpoly == charpoly2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have a characteristic polynomial, we can determine stability criteria using the `routh` function from `tbcontrol.symbolic`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tbcontrol.symbolic import routh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = routh(charpoly)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAACcAAAASCAYAAADYFMcrAAAABHNCSVQICAgIfAhkiAAAATxJREFUSInt1b0uREEYxvHfKonvRBQ0NtYlENXSaBVuQ+cClCohGhWNC6Ck0Ai1Divx0RAJhc+EAsWOhGOPXXNWQbzJ5MnM8/4zz0lO3uEX1QQWsIUbvGClCtODJZzhESeYQ3udGbsh0C32agiXx0XoW8UMNsN+H511YsAI+pFDsYZw66FnMnE+G84X68R8qmrh+oJ/jIaE14w73KMplkk2fKdGg27gOeHdYhuNGIplsoQbCFpK8Q+DFmKZLOFag16n+G/nbbFMlnDVKhf0JZbJEu7tK1tT/JZE37eZLOEOghZS/P6g7/+vGKZiFX09SvKqj4UHH0dJDBMVjh8ewrlEw3hY0I0xHCm/tXCJqXf9eeygC2vKT96g8ktTwjCuEnfEMGA6pE9bJxWYXizjHE84xTw6Kl2Qgfmvv1Gv6mF/QabV6YEAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$100$$"
      ],
      "text/plain": [
       "100"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R[0, 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All the remaining elements of the left hand row must be positive (the same sign as the first element)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "requirements = True\n",
    "for r in R[1:, 0]:\n",
    "    requirements = sympy.And(requirements, r>0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The graph below is supposed to match the textbook, but as of 2019-03-30 it does not. This appears to be a bug in `plot_implicit`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD5CAYAAAA+0W6bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHjNJREFUeJzt3XuYVNWd7vHvbprmriiCGBpvGLDlqrQCOiJ3VIwIKsooGmNEiTJmcrwk48wxkSDJkcFEjYMcJBoNiEQyYIvg4abGoIgGlYsKQoOtEJsAggrYl3X+WNXTF7qr67Kr1t5V7+d59kPVrtp7/Z62fXvXqrXX8owxiIhI+OS4LkBERBKjABcRCSkFuIhISCnARURCSgEuIhJSCnARkZBSgIuIhJQCXEQkpBTgIiIhlXUB7nneTa5rEBHxg5dtt9J7nrfTGHOy6zpERJKV67qAVPA87/2GXgJOTGctIiKpkpEBjg3pkcC+Ovs94K/pL0dExH+ZGuBFQGtjzPq6L3ietzr95YiI+C/r+sD94nleE2Ad8Jkx5jLX9YhI9snoUSie5/X3PK9NjedtPM/r59Pp7wQ2+3QuEZG4ZXSAA/8FfFXj+deRfUnxPC8fGAXMTvZcIiKJSqoL5frrCXT/y5Ilfbj00trd4C+91ItRoxoapBKb1167iu7df0Z5+UE2bZrO4MFFR71ny5ZZbN06C4CKikNcdtnGqOds1QqeeKL6+b//OxQXJ1Vmo2bPhubNU9uGiBzF8+1EyQS45wU7wGEsMAiYFHn+OLAK+O8kzlkELImcazUwPbIvmkJsd3nDjjsO9u6tft63L7z7bhJlxuDgQWjdOrVtiMhRfAvwDO9CmYkdNdgJyAfeAmYlec43gMXAqcC1wErg+iTPKSISvwwP8A7Ac8AXwN+BuZF9VaYlcM5pQAlQHDn3EODZpKoUEUlEhgd4Yxa4LsCpe+91XYGIJCPLAzzZLvxBNN7/HZt9++DOO305Vcx27kxveyLirywPcN++S/BFebnrCkQkTLI8wAM+iEZEJIoMDfCSKK+9WOPx1akuREQkZTI0wIdiR4nUNQf4cY3n/5aWaoLqzTdh5UrXVYhIojI0wB8GhgNbauybFtn/qpOKgmjPHti1y3UVIpKoDJ1O9lKgGXAJ9q7L2cDbwGvAcQ7ril1+furvxBSRcMvQK3Cw3ShPYYf6bQNWEJbwBnjoIdcViEjQZegVeBvsEEEDHMGGd4fIcw844K40ERGfZOgV+EFsSB8EvsXOInugxn6pMnMmVFa6rkJEEpGhAS6x+stfXFcgIolSgIuIhJQCXEQkpBTgwurVrisQkUQowIWZM11XICKJUICLiISUAjxANmzQre0iEjsFeIC89hps3WofN20Kbdump9316+GDD9LTloj4RwEeUKedBv/6r+lpa8uW6j8cIhIeCnARkZBSgIuIhJQCXACYOtV1BSISLwW4AHYEjIiEiwJcRCSkFOACQEUFrF3rugoRiYcCXAAoL4fHH3ddhYjEQwEeYJ06QbNmrqsQkaBSgAfYzTfbG3pEROqjAI/bYeA8oDfQHbjfbTk+mjsXVqxwXYWIxEoBHrdmwErgPWA9sBR407ezb97s26niVlZm+8JFJBwU4HHzgNaRx2WRzfPt7A884NupRCTDKcATUgH0AToAw4F+bsvx0YQJrisQkVgpwBPSBNt9UgKsBeq7jXEWUBjZStNXWpLKylxXICKxUoAnpS0wCNsPXtdEYF1ka5/GmpJz+DC8+qrrKkQkFgrwuJUC+yOPDwHLgTNT1lrHjik7db0OH4b589PbpogkRgEet13AYKAXcC62D/yylLX28MMpO7WIhFyu6wLCpxfwN9dFpNT69VBSAvn5risRkWh0BS5HWbMGtm93XYWINEYBLvXSzIQiwacAl3rNmOG6AhFpjAI8YPbu1SgQEYmNAjxgDh2Cdeuqn7dtC2embpRig0pL4dFH09+uiMROAR5wp54Ko0env92yMtixI/3tikjsFODSoKIi2LbNdRUi0hAFuDToo49sn7yIBJMCXKKaN891BSLSEAW4RLV4sesKRKQhCvAQOO88aNXKTds7dsBjj7lpW0SiU4AH0P79cORI9fOxY6FdOze1lJXBF1+4aVtEolOAB9Ds2bChvjUiHJk3D3budF2FiNSlAJdGbd0KBw64rkJE6lKAS0yuu851BSJSlwI8JPr2ddv+l1/Cxx+7rUFEalOAh8Tdd7ttf8cOeO45qKx0W4eIVFOAS8zuv19fZooEiQJc4rJmjesKRKSKAjygbrvNdQX1mzGj9hh1EXFHAR5QX35Z+/kZZ8Dw4W5qqWndOpg503UVIgIK8NBo397Nwg71efttOHzYdRUiogCXuP3xj/Dpp66rEBEFuCTk1ltdVyAiCvAQGTAAmjVzXYW1di288ILrKkSymwI8oLZsgQceqL1v/Hho08ZNPXV9/TVcey189ZXrSkSylwI8wMrLXVcQXXk5PPSQ6ypEspcCPG6fAoOBAqA78Fu35Tj2wguaI0XEFQV43HKB/wQ2A28CvwM2pa31uXPT1lRMNm6EyZM1R4qICwrwuJ0EnBN53AZ7Jf5Z2lrPCeB/sVdegQcfdF2FSPYJYByESTHwN6BfSs7+yCN2MYUw+NOfwlOrSKZQgCfsK+BK4DfAMfW8PgsojGylCbXw5Zd2TcqaunWD889P6HQp9d57cOedrqsQyS4K8ISUYcP7OmBsA++ZCKyLbO19azk/HwoKfDudr5YsgV/9ynUVItlDAR43A9yM7fv+ieNagueBB+Dll11XIZIdFOBxewN4BlgJ9IlsS1LW2vbtR+/7/vchLy9lTSbl0CEYN04LP4ikg2eMSfxgj8QPziqF2K6U+PXuDevXH72/dWt7N2RQDRtmJ73q0MF1JSKB4/l1Il2BS0osXw5TpgT7j4xI2CnAQ+qOO1xX0LjHHoMf/lDzpYikigI84PbuhQ8+OHr/BRekv5ZEPPccPPGE6ypEMpMCPOA+/RReeuno/c2aQfPm6a8nEXfdBffcc/SYdhFJjgI8pEaMgDFjXFcRu4ceggkTYMcO15WIZI5c1wVI9pg/366l+cwzwZnXXCTMNIwwLRIfRgh2Aqtt2+CUU2rvLy+Hpk2Tq8yF/Hw7xHDgQNeViDihYYTZpKGpWnNzwfPtVyF9Skpg0iR49VXXlYiEmwI8JBr6oDRjRnrr8MumTXDxxXZx5CNHXFcjEk4K8JAYP77+/SedlN46/HT4MMyaZWdYXLSo4T9SIlI/BXjIXXRRMKeXjceOHXDFFbZbZdcu19WIhIcCPCQqK+tf5LhjR+jUKf31pMITT8DZZ8OLL0JFhetqRIJPAR4Sa9fC44/X/9of/pDeWlLp73+Hyy+HsWPtl50i0jAFeAbIyYG+fV1X4a/Fi+Hcc3Ubvkg0CvAQWbas/tn98vLg9tvTX0+q7d5t+8U7doT/+A/7paeIVFOAh8iSJfDNN/W/duGF0L17eutJB2Nst8ovf2mvyJcuhc8/d12VSDAowEPm9dfr33/GGXDNNemtJd02bIBLLoF+/eB3v4PNm11XJOKWbqVPi+Rupa9pyBBYsaL+1/buhXbtfGkmFJo3hx/9CEaNgnPOgbZtXVckEhPdSp+t9u+3fcP1adUKfvCD9Nbj0uHD9k7UoUOhoABmz9ZVuWQXBXjIvPsuFBXV/1qzZnDdddm5DuXu3XDLLdC/v70af/JJDUOUzKculLTwrwsF7FzgixfbwK7P6NH29WyXl2cXhR41yn4y6dCh4Z+ZSBr51oWiAE8LfwMc4OBBuzJ9fUpLYcAA+OQTX5sMvYEDYfhwu/Xr57oayWIK8HDxP8AnTox+k8v06XD33b42mTHy8uD00+GGG+DYY+H734eWLV1XJVlEAR4u/gf44MF2Br+GVraprITLLoOXX/a12Yzkefbn2aOHvTo/80x781BDn3BEkqQADxf/Axzsqjb//M8Nv75lC3Tt6nuzWaFPHxvi3/senHqq7XLJpiGaklIK8HBJTYC3b29HX+Q0MJaorAzuvRceftj3prNOx47QogWMGwcnnmj/cFYFeq5WlpX4KMDDJTUBDlBcfPRamXUNHw7Ll6ek+azXvj1ceaWd0nfMGLuvUyfdVCRRKcDDJXUBPm6cXe09mu3b7dDDrVtTUoLUUVBgr9hzcmpPMnbhhXDCCe7qksBQgLv1A6AI6ABsiOH9qQvw5s1h3jy7ok00L78M114LBw6kpAyJQbt2dhz6McfAjTdW768aDQN2hEzTpm7qk7RRgLv1GtAauAHXAQ42FGKZavWRR2DKFNizJ2WlSJIGDLCjYAC6dLGTd1U54QQ4+WQ3dYmvFODuFQOXEYQABzvd6n33Nf6+p56yK8F/+21Ky5EU6NDBzjpZpbDQBn5No0fbL1sl0BTg7hUTPcBnRTaAUmBHSqvJy4MFC+xyZI0pKrJfvCnEM89xx9UeldSlC1x6ae33jBkD+fm197VooeBPIwW4e8UE6Qoc7P+8e/fG9t4337Qr2ivEBaBXL3sjU13DhlV36dR01lnV/fYSNwW4e8UELcBzcmDqVPjpT2N7/wsvwF132aGIIvH47nftl7F1dejQ8MIinme/sBUFeAAUE7QAB3tTyaRJ9gvLWBw6ZMct79uX2rpEIPr0BMcfb6dDjubii+0fj8a0bx/oG6wU4G6NB1YDe4ATgV8AN0d5f/oCHOzQwpkzaw9Vi+bTT+2iwU8/ndq6RNJl5MjY5rLp3t3OgxOL44+3XU0+UICHS3oDHOyY4/vvh8mTY3t/ZSX86U/wk5/AZ5+ltjaRMGrVCk47Lf7jLrjALjJSZeJEBXjIpD/AwX7JNHMmjB1rR6nE4uBBuOMOeO45fcEpkgrGaE1MicGXX8L48bBwYezHtGlju1IWL46tr1FE3FGAZ4EbbrD94fHcRj9ypB1q+MADcN55qatNRBKnLpS0cNOFUlfv3rBqlR0vHq/f/hYef9wu01ZR4X9tItlCXSiSkPfes6u2r1gR/7F33gkffQRz59pFDjThkoh7ugJPi2BcgVdp1szOTHjZZfa26iZN4j/HmjX2RqCZM+GbbyCJXyORrOLnFbgCPC2CFeA1DRoEv/+9XTYsUStXwtKlsGQJbNzoV2UimUkBHjrBDXCw05TedRcMGQLnnpv4eSor4fXX7WLLq1bB+vX+1SiSKRTgoRPsAK/SvDlMmGC7VgYPbnjF+1h88w18/rm9up8/H774Ar76Sl0tIgrw0AlHgNfUvDk8+CAMHerb7cMsWmTHpj/1FGzYAKWl/pxXJEwU4KETvgCv0qwZ9OljJxkaPLj+KUcTUVwMu3fb/vOqfvO5c/05t0iQKcBDJ7wBXlOzZnYuiKuvtgsCjBlj+889n34d9+2zY8wff9w+/+ADWL269msiYacAD53MCPD6jBljx4TffrsdjnjBBalp55VXbIjv2WMXcQbYtg127UpNeyKpogAPncwN8Lr69LFX6hMnQsuWdrx5qpSUVPejFxfDiy9Wv7ZlC/zlL6lrWyRRCvDQyZ4Ar6tFCzsrYtVKLAMG2MV4wQZ9+/apabe8vHo2xa++sqNhavrrX+Gdd2rv0zS6kg4K8NDJ3gCPpl07+Kd/so+HD69ee7FjRzvRfroVFdngr7J5MyxffvT7tm+3m0giFOChowCPxzHH2GXeqvTpAxdeWP28X7/aE+SnW2lp/UMga46oqWvjRnuTk4gCPHQU4H7LqTMN2623Hr1oxeWXwymn1N7Xtq298k+3ysqGb2LavRsWLIh+/JEjMGtWbG2VltqFOSSYFOChowAPitNPh65d63+tvj8CVfr1cxP8idiwwa5zGqu5c+Ef/4ivjSNH7CcOiZ8CPHQU4GF30kn2C9loYr3R6eqr4cQTY287iKurl5fDjh3JnePrr2HOHH/qqTJvXvx/jOIV7dNULBTgoaMAl8R4nv1kkKgf/Six6YJratMGOndO7hyZZNs2WLYs8eMnTQpIgI8Zg/n4Y9i0ya9yMpUCXMLrhBP8mw+nPlddZUcepVturl2cxIFgBDhgDhyoXmvxmWdg/357d9wLL9h9R47oFmgFuEgw1RztlC4lJQEK8Mbe8Prr9q64pUvtv5WV8P77yTQZRgpwEbEC0wdODAFeV0WFXY7r0CF48kn7RUZRUTIlhIECXESsUAd4XRUVsHevnajo2WftGNb586GszN4CnRkU4CJiZVSAN6S01M5At3OnXQhg7dowr+aiABcRKysCvK533rH95489Zm8gKClJV8v1WQrcCVQAPwR+2sj7FeAiYmVlgNf0zTf2i9CiIpg92w7crzkJUWpVAF2B/wfkA+cC84CzohyjABcRK+sDvK5Vq2x3y/LlsC7lObkG+DlQNZJ/WuTfn0U5RgEuIlZgArx79+6mRWP3F6eRMbBrVymVle05cMCOdPHfPuAAUDVL0j+Ar4GT67yvFNgTeXwE6JOKYpJQCqRoMu6EBbEmCGZdqik2QazpnY3GGF9Wl01qloUWLVqwLvWXvHEpLCxk3bp1lJfbCX3+5V/grbf8XAF9Afbqe3bk+TPAWuDRKMe0InhX4EH8VBDEmiCYdamm2ASxJu+wX2fKafwt4ZSbaxfgffFFOxfz4sXQt68fZ84Hak71VgJ8x48Ti4jEJWMDvKb27e2cB6tXw4cfwogRyZztXGALsB34FngOuDz5IkVE4pRUgE+cONGvOnwTrabWraFbNzuT2Isvwm23NT5F6NFygceAkUABMA5obP2vE+JtJA2C998umDVBMOtSTbEJYk3EuDRH45L6EnP5csywYX6V4sYbb9jpOhtaCssfQeyHExEXAjMK5dNPMV98Yec2Ads1sWyZnfx+9Gh7401ubvLzEafagQMwYwZMm1a9krm/FOAiYgUmwIkyDvyNN+xV7cCBdpz20KF2alnPi23VEleefBJ+9jM/R62AAlxEqvgZ4BhjEt7uuusu061bN9OzZ09zxRVXmH379ploVq82ZuVKYx55xJhHHzVmwQJjioqiHpKQ559/3px11lnG8zzz9ttvx338/PnGtG5tjB1Znsz2soGuBpoZmObD+fzYbjLQ3kD3ANRSte00MMjAmQbOMvCbANR0yMC5BnpFavrfAaipais30MfAqADUUrWdYqCHgd4G+gagHmNgn4ErDXQz9nfrr47r+TDy82F9ZDsA/NiYxDM44QONMSxbtsyUlZUZY4y55557zD333BNzSFZUGPP558Zs3WoD/aOPbKiXlhpz+HDcmVvLpk2bzIcffmguuuiihALcGGM2bjTmoouS+Y9VbuB0A58YOMfYINgYgF/qVw28Y4IV4J9HajIGDhj4bgB+VpUGDkYef2vgPANrAvCzMgb+08B4E7wALw1AHTW3Gwz838jjI8YGuuuajDHGADQBdgOnGJN4Bic1CmXEiBHkRlZc7d+/PyVxzDCVk2P7yrt0gcmT7Urhd9xhu1i2bIHf/96OFPn44/jrKigooFu3bvEfWMNZZ9lhh488Yochxm8tcAZwOnYFpWuBRUnV5I+BwPGui6jjJOCcyOM22NE9n7krB7D/zVpHHpdFNv8++SauBHgJO4maNOwA8Bpwc+R5HtDWXTlHGwp8YoxJamlo38aBz5kzh0suuSTp87RrZ/vIb7rJ/ltRAStWwPbtUFycfJ3xmjwZHnoI8vLiPfIzoOZKsPm4D6UwKAb+BvRzXAfYicv6AB2A4QSjph8D/4fg3cLhASOAvvg4Si4J27C30N8EnI39g/e104rquBY7C15SGr2V3vO85UB9S47eF/kowNSpU8nNzeW6665Ltp5aTjvN/ltQYEeHGAOPPmpHtjz//DBKS3cfdczUqVMZPXq0r3XceCOcfbb9d/36WI8y9ewLwhVckH0FXAn8BjjGcS1gP+WuB/YDY4ANgMtv4Iuwf0z6Aqsd1lGfN7B3JH+B/WN3JvbTnivlwLvYKS76Yad//hUwxWFNlud5edi7/6LNgBeTRgPcGBN1pPfTTz9NUVERK1aswPNSF1BVV8CTJ9upY885Zzn79sHFF6esyVp69YInnoCbb4YNG2I5Qrfcx6cMG97XAWMd11JXW2AQdh54lwH+BrAYWAIcxnYTXA8867CmKlW/2x2wf+zW4jbA8yNb1aemq7ABHgiXAO8aY/6e7ImS+hy2dOlSfv3rX7N48WJatmyZbC0xy82Ffv2qw3vp0vS0e9558Ic/VH8yiK7mLfcG3XIfjcH2VRYAP3FcS5VS7JU3wCFgOfaq0qVp2AuBYuzv0xCCEd5fAwdrPH4Ft3/owHYadAY+ijxfQfQ5+9NqPD50nwAk/O2nMYYuXbqY/Px807t3b9O7d29z6623JjTiwy/z5hnz1lvGLFy40HTq1Mnk5eWZDh06mBEjRvjazqJFxhx3XCzfOL9k7IiKPAO/dP7tt92uNdDRQK6BTgZmB6Cm1w1goKexw6x6R352Lmt6z9ihej2NHbHziwD8nGpuq0xwRqF8Yuwoq6ohl0H5Xf+bsUMaexoYbWBvAGr62mDnoD7WmMSzt2rLiAUdXFi5EsaNs6sBNU438oiIZXy8kSdoX2X76umn4YMPUnPuIUPg5z9PzblFRGKR0VfgVR9aclL0Z6q8HK6/HubPb+ydugIXEUtX4DHyvNrhHfsQwNjk5tox4n2CtlqaiGSFjA7wujrWN5o9SZ07w7NBGAggIllHAe6DLl3sF5oiIumUVQGeKs2bw/Tp0L2xhXlERHykAPdJ585www2uqxARd1rXeLwE+C6wM+ajPc/r7HneKs/zNnuet9HzvDsbO0YB7qPbb4errnJdhYi4tQKYjJ164eR4DiwH/pcxpgDoD9zueV7U20cV4D5q1QomTIBjj3VdiYi48TpwC3bK3y4NvsvzvBM9z/uz53nvRbbzjTG7jDHvAhhjDgKbgU7RWlOA++zyy6FnT9dViEj6HQFGA/9NDPPmPAK8aozpjZ0Mv9ay6p7nnYqdB/etaCdRgKfAk0+6rkBE0q8pcD4QUwAMAf4LwBhTYYz5suoFz/NaAy9gl1s7EO0kCvAUOPlkuOUW11WISHrlAM8DbwMPJnQGz/OaYsP7j8aYhbG0KD5r3twuEZeqW/hFJKhaYhfe+CONXImvACYBeJ7XxPO8Yzy7oMKTwGZjzIxYWsvouVBcKyiADz8EzYUikg1aY1eVAruYy0Ds6lK1VwgzBs/zvBOxa8+djl27bxJ2CajXgQ+Aysjb/80Ys6ShFhXgKbRoEVxxBSjARaSKJrMKiTPOgDZtXFchIplKAR6HBQsW0L17d3Jycli3rvEr6u7dYdKkNBQmIgE1FehTa/M87z6/zq4Aj0OPHj1YuHAhAwfGvljrgAH6MlMke90HrK+1GWOm+nX2Rlell2oFBQVxH3PFFXbe8G+/TUFBIpLVdG2YBqmaxlZEspuuwOsYNmwYu3fvPmr/1KlTGT16dD1H1G/WrFnMmjULgCNHSn2rT0SkioYRJmDQoEFMnz6dwsLCmN7fu3ch7dqtY9WqFBcmIoGnYYQh07Qp9OjhugoRyTQK8Dj8+c9/Jj8/nzVr1jBq1ChGjhwZ87FTpsBxx6WwOBHJOupCSYPCwkLWrVtH+/awZ4/rakTEJXWhhNQvfuG6AhHJJArwNOra1XUFIpJJFOBx2rsXShMcFdijBwwa5Gs5IpLFFOBxysmBJk0SO7ZjR+jS8DJ5IiJx0Y08cWrb1u3xIiJVdAWeZtOnQ8uWrqsQkUygAE/Q8uWuKxCRbKcAT9DgwYkfO2GCf3WISPZSgCco0S8yAUaN8q8OEcleCvAEffwxfPKJ6ypEJJspwBPUtWviQwK/9z0YP97fekQk+yjAk7BsmesKRCSbKcCTEMdkhEdp2tS/OkQkOynAHXn6adcViEjYKcCTUF6u8eAi4o4CPAm5uTBsWOLH645MEUmGAtyhOXNcVyAiYaYAT9LmzbBzp+sqRCQbaTbCJBUUuK5ARLKVrsAdOv98KCx0XYWIhJUC3Af33pvYcZ07Q6dO/tYiItlDAe6DKVNcVyAi2UgB7oO8vMSP7dnTvzpEJLsowB277z7XFYhIWCnAfXDoEKxe7boKEck2CnAftGgBgwa5rkJEso0C3LEmTWD4cNdViEgYKcDjcPfdd3PmmWfSq1cvxowZw/79+5M+Z9OmWtxBRBKjAI/D8OHD2bBhA++//z5du3Zl2rRp//Pa++/D5587LE5Eso4CPA4jRowgN9fOPtC/f39KSkr+57VeveA733FVmYhkI82FkqA5c+ZwzTXX+HKudu3gnHN8OZWIZBHPGOO6hkDxPG850LGel+4zxiyKvOc+oBAYaxr4AXqeNxGYGHna3BjTIxX1ikj2UoDHyfO8G4HbgKHGmG9c1yMi2UtdKHHwPO9i4F7gIoW3iLimK/A4eJ63FWgG/COy601jzG0OSxKRLKYAFxEJKQ0jFBEJKQW4iEhIKcBFREJKAS4iElIKcBGRkFKAi4iElAJcRCSkFOAiIiH1/wF19lNNhADmJgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<sympy.plotting.plot.Plot at 0x1217574e0>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sympy.plot_implicit(requirements, (K_c2, -2, 7), (K_c1, -2, 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an alternative, let's evaluate numerically"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = sympy.lambdify((K_c2, K_c1), requirements)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "nK_c2, nK_c1 = numpy.meshgrid(numpy.linspace(-2, 4, 300), numpy.linspace(-2, 7, 300))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = f(nK_c2, nK_c1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'K_c2')"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAELCAYAAAA2mZrgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAELNJREFUeJzt3X2MZXV9x/H3p8sCIrsxtBZZlobaKJESV9oVUFK1QN1VqU8pDbSa+tBM2mgLjYkV+ac2sf2jiaFJje1WrX9I1fpAtPiwrgqlpoouukvBRUWicV3b1agVJFkW+PaPuXsYl9mZOzv3zu/ec9+vZLJz7z1z53P26TO/3/mdc1JVSJIE8AutA0iSJoelIEnqWAqSpI6lIEnqWAqSpI6lIEnqNC2FJOck2bPg46dJrmmZSZJmWSblPIUk64DvARdW1Xda55GkWTRJ00eXAt+yECSpnRNaB1jgSuB9i72QZA6YA1jHut88hY1rmWvsnvr0B1pHmGjfuOOU1hGkqXcfP/5hVT1xue0mYvooyYnAAeDXq+p/l9p2Y06rC3Pp2gRbIzsP7G0dYaJt27SldQRp6n2mPnR7VW1dbrtJmT56AfCV5QpBkjRek1IKV3GMqSNJ0tppXgpJTgF+B/hI6yySNOuaH2iuqgeAX2ydQ5I0ASMFSdLksBQkSR1LQZLUsRQkSR1LQZLUsRQkSR1LQZLUsRQkSR1LQZLUsRQkSR1LQZLUsRQkSR1LQZLUsRQkSR1LQZLUaV4KSZ6Q5ENJ7k6yL8mzWmeSpFnV/CY7wN8Dn6qq30tyInBK60CSNKualkKSjcBzgFcBVNWDwIMtM0nSLGs9ffRk4AfAvyT5apJ3Jnl840ySNLNal8IJwG8A76iq84GfAW86eqMkc0l2J9l9mENrnVGSZkbrUtgP7K+q2waPP8R8SfycqtpRVVuraut6TlrTgJI0S5qWQlX9D/DdJOcMnroU+FrDSJI00yZh9dGfATcMVh7dC7y6cR5JmlnNS6Gq9gBbW+eQJLU/piBg26YtrSNIEmApSJIWsBQ08XYe2Ns6gjQzLAVJUsdSkCR1LAVJUsdSkCR1LAVJUsdSkCR1LAVJUsdSkCR1LAVJUsdSkCR1LAVJUsdSkCR1LAVJUqf5TXaSfBu4D3gYeKiqvOGOJDXSvBQGfruqftg6hCTNOqePJEmdSSiFAj6d5PYkc4ttkGQuye4kuw9zaI3jSdLsmITpo4ur6kCSXwZ2Jbm7qm5duEFV7QB2AGzMadUipCTNguYjhao6MPj1IHAjcEHbRJI0u5qWQpLHJ9lw5HPg+cCdLTNJ0ixrPX10OnBjkiNZ/rWqPtU2kiTNrqalUFX3AltaZpAkPar5MQVJ0uSwFCRJHUthQmzb5CyapPYsBU2FnQf2to4gzQRLQZLUsRQkSR1LQZLUsRQkSR1LQZLUsRQkSR1LQZLUsRQkSR1LQZLUsRQkSR1LQZLUmYhSSLIuyVeT3NQ6iyTNsokoBeBqYF/rEJI065qXQpLNwIuAd7bOIkmzrnkpANcDbwQeOdYGSeaS7E6y+zCH1i6ZJM2YpqWQ5HLgYFXdvtR2VbWjqrZW1db1nLRG6SRp9rQeKVwMvDjJt4H3A5ckeW/bSJI0u5qWQlVdW1Wbq+ps4Ergc1X1ipaZWvKWnEvz7mvS+LUeKUiSJsgJrQMcUVW3ALc0jiFJM82RgiSpYylIkjqWgiSpYylIkjqWgiSpYylIkjqWgiSpYylIkjqWgiSpYyloqnj9I2m8LAVJUsdSmDDbNm3xaqmSmrEUJEkdS0GS1LEUJEmd1vdoPjnJl5LsTXJXkre0zCNJs671TXYOAZdU1f1J1gOfT/LJqvpi41ySNJOalkJVFXD/4OH6wUe1SyRJs635MYUk65LsAQ4Cu6rqtkW2mUuyO8nuwxxa+5CSNCOal0JVPVxVzwA2AxckOW+RbXZU1daq2rqek9Y+pCaKZzVL49O8FI6oqp8AtwDbG0eRpJnVevXRE5M8YfD544DLgLtbZpKkWbaqUkjy6lV+/zOAm5PcAXyZ+WMKN63yPSVJx2m1I4VVnVdQVXdU1flV9fSqOq+q/nqVeXrD6x9JamHZJamDn+IXfQk4fbRxJEktDXOewunANuDHRz0f4L9GnkiS1Mww00c3AadW1XeO+vg286uFpDXnslRpPJYdKVTVa5d47Q9GG0eS1NLQB5qTXJRkw4LHG5JcOJ5YkqQWVrL66B08ep0igJ8NnpMk9cRKSiGDC9gBUFWP0P4qq5KkEVpJKdyb5M+TrB98XA3cO65gkqS1t5JS+BPg2cD3gP3AhcDcOEJpniewLW3ngb2uQpJGbOjpn6o6CFx5rNeTXFtVfzuSVJKkJkZ5QbwrRvhekqQGRlkKGeF7SZIaGGUpeBtNSZpyjhQ09TzYLI3OsqWQZPMSr/3ugocfHEkiSVIzw4wUPpvk7KOfTPIa4Pojj6vqb1b6zZOcleTmJPuS3DU490GS1MgwpfAXwK4kTznyRJJrB88/d5Xf/yHgDVX1NOAi4HVJzl3le/aK5yoMxykkaTSGuUrqJ5IcAj6Z5KXAHwPPBJ5TVUffY2FFqur7wPcHn9+XZB9wJvC11byvJOn4DHWguao+C7yK+fsnPBm4dLWFcLTBFNX5wG2LvDaXZHeS3Yc5NMpvK0laYJjbcd7H/HLTACcBlwIHkwSoqtq42hBJTgU+DFxTVT89+vWq2gHsANiY01z6KkljsuxIoao2VNXGwa8nVtXjFzweRSGsZ74Qbqiqj6z2/TS7PK4grd4oz1NYscFo413Avqp6W8sskqTGpQBcDLwSuCTJnsHHCxtnmjiuQJK0VpreJKeqPo9nQmuEdh7Ya4lKq9B6pCBJmiCWgiSpYymod1yFJB0/S2FKOE8uaS1YCuolRwvS8bEUJEkdS0G95WhBWjlLQZLUsRSmiAebV27ngb2OGKQVsBQ0EywGaTiWgiSpYyloZjhakJZnKUwZjyusjsUgLc1SkCR1LAXNHFckScfWvBSSvDvJwSR3ts4iSbOueSkA7wG2tw4xTbZt2uKxhRFwxCA9VvNSqKpbgR+1zqHZZTlIj2peCsNIMpdkd5LdhznUOo56ymKQpqQUqmpHVW2tqq3rOal1nInhFNLoOWrQrJuKUpDWmsWgWWUpSMdgMWgWNS+FJO8DvgCck2R/kte2zjRNnEIaryPTSRaEZsUJrQNU1VWtM0jDWFgMlrH6qnkpaPW2bdriT7Jr7Ojfb0tCfdF8+kjqA6eY1BeOFKQRcopJ085S6AmnkCbPsf48LAtNMktBWmOLlYVFoUlhKfSIo4Xp5YFrTQpLoWcshn44nj9Di0SjYClIPbFckVgaGoal0EOOFrSYUf+dsGT6yVKQdFyOt2Qsk8lmKfTUkX94jhg0aRyxTDZLoeecSlLf+ff75622JC0FSeqRY5XkujOG+3qvfTQDHF5LGpalMCMsBknDaF4KSbYn+XqSe5K8qXWePrMYJC2naSkkWQe8HXgBcC5wVZJzW2bqO4tB0lJajxQuAO6pqnur6kHg/cBLGmfqPYtB0rG0LoUzge8ueLx/8JzGzGKQtJjWS1KzyHP1mI2SOWAO4GROGXemmbGwGFzrLQnajxT2A2cteLwZOHD0RlW1o6q2VtXW9Zy0ZuEkada0LoUvA09J8qtJTgSuBD7WONNMcjpJEjSePqqqh5K8HtgJrAPeXVV3tcw0y5xOktR6pEBVfaKqnlpVv1ZVb22dR/McOUizqXkpaHJt27TFcpBmTOvVR5oCTitJs8ORglbEkYPUb44UtGJHF4OjB6k/HClo1Tz2IPWHIwWNjMcepOlnKWgsFhs5WBTS5LMUtGYcSUiTz1JQE0sdg7AwpHYsBU0cp56kdiwFTQWXwUprw1LQVFrpElhLRBqOpaCZcKwSsSykn2cpaKYtN+KwNDRrLAVpCeM4U9ui0SRrVgpJrgD+CngacEFV7W6VRVpLx1s0lonWQsuRwp3Ay4F/aphBmhqTen0py6pfmpVCVe0DSNIqgqQRWOuysoTGy2MKkqbKJI6Y+lRUYy2FJJ8BnrTIS9dV1UdX8D5zwBzAyZwyonSSNBqTWFSP9c2hthprKVTVZSN6nx3ADoCNOa1G8Z6SpMfyJjuSpE6zUkjysiT7gWcBH0+ys1UWSdK8lquPbgRubPX9JUmP5fSRJKljKUiSOpaCJKljKUiSOpaCJKljKUiSOpaCJKljKUiSOpaCJKljKUiSOpaCJKljKUiSOpaCJKljKUiSOpaCJKljKUiSOi3vvPZ3Se5OckeSG5M8oVUWSdK8liOFXcB5VfV04BvAtQ2zSJJoWApV9emqemjw8IvA5lZZJEnzUlWtM5Dk34EPVNV7j/H6HDA3eHgecOdaZWvgl4Aftg4xJn3eN3D/pl3f9++cqtqw3EZjLYUknwGetMhL11XVRwfbXAdsBV5eQ4RJsruqto426eTo8/71ed/A/Zt27t+8E8YZoqouW+r1JH8EXA5cOkwhSJLGa6ylsJQk24G/BJ5bVQ+0yiFJelTL1Uf/AGwAdiXZk+Qfh/y6HWPMNAn6vH993jdw/6ad+8eEHGiWJE0Gz2iWJHUsBUlSZypLoc+XyEhyRZK7kjySpDfL45JsT/L1JPckeVPrPKOU5N1JDibp5fkzSc5KcnOSfYO/m1e3zjQqSU5O8qUkewf79pbWmcYhybokX01y03LbTmUp0O9LZNwJvBy4tXWQUUmyDng78ALgXOCqJOe2TTVS7wG2tw4xRg8Bb6iqpwEXAa/r0Z/fIeCSqtoCPAPYnuSixpnG4Wpg3zAbTmUp9PkSGVW1r6q+3jrHiF0A3FNV91bVg8D7gZc0zjQyVXUr8KPWOcalqr5fVV8ZfH4f8/+5nNk21WjUvPsHD9cPPnq1+ibJZuBFwDuH2X4qS+EorwE+2TqElnQm8N0Fj/fTk/9UZk2Ss4HzgdvaJhmdwdTKHuAgsKuqerNvA9cDbwQeGWbjZievLWcFl8h4CLhhLbOt1jD71jNZ5Lle/TQ2C5KcCnwYuKaqfto6z6hU1cPAMwbHJm9Mcl5V9eL4UJLLgYNVdXuS5w3zNRNbCn2+RMZy+9ZD+4GzFjzeDBxolEXHIcl65gvhhqr6SOs841BVP0lyC/PHh3pRCsDFwIuTvBA4GdiY5L1V9YpjfcFUTh8tuETGi71ExlT4MvCUJL+a5ETgSuBjjTNpSEkCvAvYV1Vva51nlJI88cjqxSSPAy4D7m6banSq6tqq2lxVZzP/7+5zSxUCTGkpcPyXyJh4SV6WZD/wLODjSXa2zrRag0UBrwd2Mn+Q8t+q6q62qUYnyfuALwDnJNmf5LWtM43YxcArgUsG/972DH7y7IMzgJuT3MH8Dy+7qmrZZZt95mUuJEmdaR0pSJLGwFKQJHUsBUlSx1KQJHUsBUlSx1KQJHUsBWkZSe5f8PkLk3wzya+s4Ot7e+lp9Y/nKUjLSHJ/VZ2a5FLm73P7/Kr61gq+/gzgjKr6SpINwO3AS6vqa2OKLB03RwrSEJL8FvDPwIuWKoQkpw9u/LR38PHsPl96Wv3jSEFaRpLDwH3A86rqjmW2/QDwhaq6fnBzoVOr6v8WvH428zdQOq9PVxpVf1gK0jKSPAB8DvhWVS15PCDJD4DNVXVokddOBf4DeGtfrzSq6ef0kbS8R4DfB56Z5M3H8wazcOlp9YOlIA1hcIn2y4E/XOYqqJ8F/hS6O3pt7POlp9U/Th9Jyziy+mjw+VnMHxO4ZrG75CU5nfkVSk8GHma+INYB/wn8N4/eEvHNVfWJNYgvrYilIEnqOH0kSepM7D2apUmW5DrgiqOe/mBVvbVFHmlUnD6SJHWcPpIkdSwFSVLHUpAkdSwFSVLn/wFF4u7OaNxtEAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.pcolor(nK_c2, nK_c1, r)\n",
    "plt.ylabel('K_c1')\n",
    "plt.xlabel('K_c2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that even this simple system can exhibit more complicated behaviour than we may expect from first order systems because of the extra loops formed by the controllers."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}